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GRADIENT SCHEMES FOR THE SIGNORINI AND THE 
OBSTACLE PROBLEMS, AND APPLICATION TO HYBRID 
MIMETIC MIXED METHODS 

YAHYA ALNASHRI AND JEROME DRONIOU 


Abstract. Gradient schemes is a framework which enables the unified conver¬ 
gence analysis of many different methods - such as finite elements (conforming, 
non-conforming and mixed) and finite volumes methods - for 2 nd order dif¬ 
fusion equations. We show in this work that the gradient schemes framework 
can be extended to variational inequalities involving mixed Dirichlet, Neu¬ 
mann and Signorini boundary conditions. This extension allows us to provide 
error estimates for numerical approximations of such models, recovering known 
convergence rates for some methods, and establishing new convergence rates 
for schemes not previously studied for variational inequalities. The general 
framework we develop also enables us to design a new numerical method for 
the obstacle and Signorini problems, based on hybrid mimetic mixed schemes. 
We provide numerical results that demonstrate the accuracy of these schemes, 
and confirm our theoretical rates of convergence. 


1. Introduction 

Numerous problems arising in fluid dynamics, elasticity, biomathematics, mathe¬ 
matical economics and control theory are modelled by second order partial differ¬ 
ential equations with one-sided constraints on the solution or its normal derivative 
on the boundary. Such models also appear in free boundary problems, in which 
partial differential equations (PDEs) are written on a domain whose boundary is 
not given but has to be located as a part of the solution. From the mathematical 
point of view, these models can be recast into variational inequalities 03 El 1111- 
Linear variational inequalities are used in the study of triple points in electrochem¬ 
ical reactions [35] , of contact in linear elasticity [33] , and of lubrication phenomena 
m- Other models involve non-linear variational inequalities; this is for example 
the case of unconfined seepage models, which involve non-linear quasi-variational 
inequalities obtained through the Baiocchi transform 0], or quasi-linear classical 
variation inequalities through the extension of the Darcy velocity inside the dry 
domain 051 [El [361 S3- It is worth mentioning that several of these linear and non¬ 
linear models involve possibly heterogeneous and anisotropic tensors. For linear 
equations, this happens in models of contact elasticity problems involving com¬ 
posite materials (for which the stiffness tensor depends on the position), and in 
lubrication problems (the tensor is a function of the first fundamental form of the 
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film). For non-linear models, heterogeneity and anisotropy are encountered in seep¬ 
age problems in porous media. 

In this work, we consider two types of linear elliptic variational inequalities: the 
Signorini problem and the obstacle problem. The Signorini problem is formulated 


as 



(i.i) 

—div(AVu) = / 

in 12, 

(1.2) 

u = 0 

on Ti, 

(1.3) 

AVw•n = 0 

on r 2 , 

(1.4) 

u< a ) 

AVu • n <0 > 

AVii • n (a — u) = 0 J 

on r 3 . 

The obstacle problem 

is 


(1.5) 

(div(AVu) + f){g-u) = 0 

in 0, 

(1.6) 

—div(AVw) < / 

in 12, 

(1.7) 

u<g 

in fi, 

(1.8) 

u = 0 

on 9S2 


Here H is a bounded open set of (d > 1), n is the unit outer normal to dfl and 
(ri,r 2 ,r 3 ) is a partition of dfl (precise assumptions are stated in the next section). 

Mathematical theories associated with the existence, uniqueness and stability of 
the solutions to variational inequalities have been extensively developed [21, ;34] [29]. 
From the numerical point of view, different methods have been considered to ap¬ 
proximate variational inequalities. Bardet and Tobita [5] applied a finite differ¬ 
ence scheme to the unconfined seepage problem. Extensions of the discontinuous 
Galerkin method to solve the obstacle problem can be found in mm- Although 
this method is still applicable when the functions are discontinuous along the el¬ 
ements boundaries, the exact solution has to be in the space H 2 to ensure the 
consistency and the convergence. Numerical verification methods, which aim at 
finding a set in which a solution exists, have been developed for a few variational 
inequality problems. In particular we cite the obstacle problems, the Signorini 
problem and elasto-plastic problems (see |.'I8| and references therein). 

Falk [26] was the first to provide a general error estimate for the approximation 
by conforming methods of the obstacle problem. This estimate showed that the Pi 
finite element error is 0(h). Yang and Chao [3Dj showed that the convergence rate 
of non-conforming finite elements method for the Signorini problem has order one, 
under an H 5 / 2 -regularity assumption on the solution. 

Herbin and Marchand [32] showed that if A = Id and if the grid satisfies the 
orthogonality condition required by the two-point flux approximation (TPFA) finite 
volume method, then for both problems the solutions provided by this scheme 
converge in L 2 (fl) to the unique solution as the mesh size tends to zero. 

The gradient schemes framework is a discretisation of weak variational formu¬ 
lations of PDEs using a small number of discrete elements and properties. Previ¬ 
ous studies [19] showed that this framework includes many well-known numerical 
schemes: conforming, non-conforming and mixed finite elements methods (including 
the non-conforming “Crouzeix- Raviart” method and the Raviart-Thomas method), 
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hybrid mimetic mixed methods (which contain hybrid mimetic finite differences, hy¬ 
brid finite volumes/SUSHI scheme and mixed finite volumes), nodal mimetic finite 
differences, and finite volumes methods (such as some multi-points flux approxi¬ 
mation and discrete duality finite volume methods). Boundary conditions can be 
considered in various ways within the gradient schemes framework, which is a use¬ 
ful flexibility for some problems and engineering constraints (e.g. available meshes, 
etc.). This framework has been analysed for several linear and non-linear elliptic 
and parabolic problems, including the Leary-Lions, Stokes, Richards, Stefan’s and 
elasticity equations. It is noted that, for models based on linear operators as in 
(1.1) only three properties are sufficient to obtain the convergence of schemes for 
the corresponding PDEs. We refer the reader to for more 

details. 

The aim of this paper is to develop a gradient schemes framework for elliptic lin¬ 
ear variational inequalities with different types of boundary conditions. Although 
we focus here on scalar models, as opposed to variational inequalities involving dis¬ 
placements, we notice that using [20] our results can be easily adapted to elasticity 
models. The gradient scheme framework enables us to design a unified convergence 
analysis of many numerical methods for variational inequalities; using the theo¬ 
rems stemming from this analysis, we recover known estimates for some schemes, 
and we establish new estimates for methods which were not, to our best knowl¬ 
edge, previously studied for variational inequalities. In particular, we apply our 
results to the hybrid mimetic mixed (HMM) method of 16], which contains the 
mixed/hybrid mimetic finite difference methods and, on certain grids (such as tri¬ 
angular grids), the TPFA scheme. We would like to point out that, to our best 
knowledge, although several papers studied the obstacle problem for nodal mimetic 
finite differences (see e.g. [ 2113 ]), no study has been conducted of the mixed/hybrid 
mimetic method for the obstacle problem, or of any mimetic scheme for the Sig- 
norini problem. Our analysis through the gradient schemes framework, and our 
application to the HMM method therefore seems to provide entirely new results for 
these schemes on meaningful variational inequality models. We also mention that 
we are mostly interested in low-order methods, since the models we consider could 
involve discontinuous data - as in heterogeneous media - which prevent the solution 
from being smooth. Finally, it seems to us that our unified analysis gives simpler 
proofs than the previous analyses in the literature. The study of linear models is a 
first step towards a complete analysis of non-linear models, a work already initiated 

in D-. 

This paper is organised as follows. In Section [2] we present the gradient schemes 
framework for variational inequalities, and we state our main error estimates. We 
then show in Section [3] that several methods, of practical relevance for anisotropic 
heterogeneous elliptic operators, fit into this framework and therefore satisfy the 
error estimates established in Section [2] This also enables us to design and analyse 
for the first time an HMM scheme for variational inequalities. We provide in Section 
|4] numerical results that demonstrate the excellent behaviour of this new scheme 
on test cases from the literature. Available articles on numerical methods for vari¬ 
ational inequalities do not seem to present any test case with analytic solutions, on 
which theoretical rates of convergence can be rigorously checked; we therefore de¬ 
velop in Section [4] such a test case, and we use it to assess the practical convergence 
rates of HMM for the Signorini problem. Section [5] is devoted to the proof of the 
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main results. Since we deal here with possibly nonconforming schemes, the tech¬ 
nique used in [25] cannot be used to obtain error estimates. Instead, we develop a 
similar technique as in [23] for gradient schemes for PDEs. Dealing with variational 
inequalities however requires us to establish new estimates with additional terms. 
These additional terms apparently lead to a degraded rate of convergence, but we 
show that the optimal rates can be easily recovered for many classical methods. We 
present in Section [ 6 ] an extension of our results to the case of approximate barriers, 
which is natural in many schemes in which the exact barrier is replaced with an 
interpolant depending on the approximation space. The paper is completed by two 
short sections: a conclusion, and an appendix recalling the definition of the normal 
trace of H^w functions. 


2. Assumptions and main results 

2.1. Weak formulations. Let us start by stating our assumptions and weak for¬ 
mulations for the Signorini and the obstacle problems. 

Hypothesis 2.1 (Signorini problem). We make the following assumptions on the 
data in 

(1) 12 is an open bounded connected subset ofM. d , d> 1 and fI has a Lipschitz 
boundary, 

(2) A is a measurable function from 12 to M^(K) (where M^R) is the set of 
dx d matrices) and there exists X, A > 0 such that, for a.e. x £ 12, A(x) is 
symmetric with eigenvalues in [A, A], 

(3) Ti,T 2 andT 3 are measurable pairwise disjoint subsets of dQ such that TiL) 
T 2 U T 3 = 912 and Ti has a positive (d — 1)-dimensional measure, 

(4) /eL 2 ( 12 ), aeL 2 (912). 


Under Hypothesis 2.1 the weak formulation of Problem (1.1)- (1.4) is 

Find u £ K. := {v £ 4L 1 (D) : 7(4;) = 0 on 7 ( 1 ;) < a on r 3 } s.t., 

' Vt £ K ,, f A(x)Vu(x) ■ V(m — v)(x) dx < f f(x)(u(x) — v(x))dx, 
Jn J n 


where 7 : H 1 ( 12) ha it 1 / 2 (912) is the trace operator. We refer the reader to (2S] 
for the proof of equivalence between the strong and the weak formulations. It was 
shown in [ 21 ] that, if K, is not empty (which is the case, for example, if a > 0 on 
912) then there exists a unique solution to Problem & In the sequel, we will 
assume that the barrier a is such that K, is not empty. 


Hypothesis 2.2 (Obstacle problem). 
are: 


Our assumptions on the data in (1.5) (1.8) 


(1) 12 and A satisfy (1) and (2) in Hypothesis 2.1 

(2) / £ L 2 (Q), g £ L 2 (fl). 


Under Hypotheses 2.2 we consider the obstacle problem (1.5)- (1.8) under the 
following weak form: 

Find u £ /C := {v £ Hg(il) : v < g in 12} such that, for all v £ /C, 

/ A(x)'\7u(x) ■ V(m — v)(x) da; < / f(x)(u(x) — v(x))dx. 

Jn Jn 


( 2 . 2 ) 











GRADIENT SCHEMES FOR THE SIGNORINI AND THE OBSTACLE PROBLEMS 


5 


It can be seen [43] that if u £ C 2 (d) and A is Lipschitz continuous, then (1.5)- (1.8) 
and (2.2) are indeed equivalent. The proof of this equivalence can be easily adapted 
to the case where the solution belongs to H 2 (Cl). If g is such that 1 C is not empty, 
which we assume from here on, then ( 2 . 2 ) has a unique solution 


2.2. Construction of gradient schemes. Gradient schemes provide a general 
formulation of different numerical methods. Each gradient scheme is based on 
a gradient discretisation, which is a set of discrete space and operators used to 
discretise the weak formulation of the problem under study. Actually, a gradient 
scheme consists of replacing the continuous space and operators used in the weak 
formulations by the discrete counterparts provided by a gradient discretisation. 
In this part, we define gradient discretisations for the Signorini and the obstacle 
problems, and we list the properties that are required of a gradient discretisation 
to give rise to a converging gradient scheme. 


2.2.1. Signorini problem. 

Definition 2.3. (Gradient discretisation for the Signorini problem). A gradient 
discretisation T> for Problem (2.1) is defined by V = (1^^,lip,Tp, Vp), where 
(1) the set of discrete unknowns X-p l r 1 is a finite dimensional vector space on 
taking into account the homogeneous boundary conditions on Pi, 
the linear mapping lip : Xx>Xi —► L 2 (Cl) is the reconstructed function, 
the linear mapping Tp : Ipy, —> H v ' (<9f2) is the reconstructed trace, 
where i/jV 2 (<9fl) = {u € H 1 ^ 2 (dCl) : v = 0 on Pi}, 

the linear mapping Vp : Ip^i —► L 2 (Ci) d is a reconstructed gradient, 
which must be defined such that ||Vp • is a norm on Xp^ 1 . 

As explained above, a gradient scheme for 413 is obtained by simply replacing 
in the weak formulation of the problem the continuous space and operators by the 
discrete space and operators coming from a gradient discretisation. 

Definition 2.4. (Gradient scheme for the Signorini problem). Let V be the gra¬ 
dient dicretisation in the sense of Definition |2.3| The gradient scheme scheme for 
Problem (2.1) is 


( 2 ) 

(3) 

(4) 


(2.3) 


Find u £ AC p := {v £ Xp^ x '■ Tpv < a on Ta} such that, \/v £ AC p, 

A(i)Vpi/(a;) • Vp(it — v)(x) da; < / f(x)Hp(u — v)(x)dx. 

’ J n 


The accuracy of a gradient discretisation is measured through three indicators, 
related to its coercivity, consistency and limit-conformity. The good behaviour 
of these indicators along a sequence of gradient discretisations ensures that the 
solutions to the corresponding gradient schemes converge towards the solution to 
the continuous problem. 

To measure the coercivity of a gradient discretisation T> in the sense of Definition 
2.3[ we define the norm Cp of the linear mapping lip by 

( ||npV||p2(Q) i ||Tp1l||p-l/2( a Q) \ 

l|Vpu|| iW / 


(2.4) 


Cp = 


max 


#6^o, ri \{0} V||Vp^|U2 (0)d 
The consistency of a gradient discretisation V in the sense of Definition |2.3| is 
measured by <Sp : K. x /Cp —► [0, +oo) defined by 

(2.5) V(vj,u) £ /Cx/Cp, Sp(ip,v) = Wllpv - ip\\ L 2( n) + \\Xpv -Vip\\ L 2( Q y. 
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To measure the limit-conformity of the gradient discretisation T> in the sense of 
Definition 2.3 we introduce Wp '■ H dd ,(Q) 


[0, +oo) defined by 


(2.6) Vi/> € i/ d iv(fi) 

Wp(t/>) = 


sup 

»EX B , ri \{0} II V-DJ’|| L 2( n) d 


(\7' D v-ip+H'DV-div(ip)) dec-( 7 „(■»/>), T-pu) 


where H div (fl) = {xf £ L 2 (tt) d : div0 £ L 2 (D)}. We refer the reader to the 
appendix (Section [s]) for the definitions of the normal trace 7 „ on H d - lv (Q), and 
of the duality product (•, •) between (K3 rl / 2 (e)Q)y and H 1 ^ 2 (d^l). 

In order for a sequence of gradient discretisations (V m ) me ^ to provide converging 
gradient schemes for PDEs, it is expected that (Cp iTi ) m6 N remains bounded and that 
(inf V £K. T > rn Sp m (•, v))me n and (Wc m )me n converge pointwise to 0. These properties 
are precisely called the coercivity, the consistency and the limit-conformity of the 
sequence of gradient discretisations (I2 m ) m eNj see [18] . 

The definition (2.4) of Cp does not only include the norm of lip, as in the 


obstacle problem case detailed below, but also quite naturally the norm of the 
other reconstruction operator Tp. 


The definition (2.6) takes into account the non-zero boundary conditions on 
a part of dfl. This was already noticed in the case of gradient discretisations 
adapted for PDEs with mixed boundary conditions, but a major difference must 
be raised here. Wp will be applied to if = AVu. For PDEs [TSJ iT5] the 
boundary conditions ensure that 7 n (i/>) £ L 2 {d£l) and thus that (7n(■*/>), Tpw) can 
be replaced with f r ^ 7„(’0)Tpj; da; in Wp. Hence, in the study of gradient schemes 
for PDEs with mixed boundary conditions, Tpi; only needs to be in L 2 (dQ). For 
the Signorini problem we cannot ensure that 7 n (i/>) £ L 2 (d H); we only know that 
7n(V0 € U 1 / 2 (9 fl). The definition of Tp therefore needed to be changed to ensure 
that this reconstructed trace takes values in H 1 ^ 2 (dfl) instead of L 2 (dfl). 

2.2.2. Obstacle problem. Gradient schemes for the obstacle problems already ap¬ 
pear in PQ. We just recall the following definitions for the sake of legibility. 

The definition of a gradient discretisation for the obstacle problem is not different 
from the definition of a gradient discretisation for elliptic PDEs with homogeneous 
Dirichlet conditions [18] [21]. 

Definition 2.5. (Gradient discretisation for the obstacle problem). A gradient 
discretisation T> for homogeneous Dirichlet boundary conditions is defined by a 
triplet D = (Xp j0 ,Dp, Vp), where 

(1) the set of discrete unknowns Ap j0 is a finite dimensional vector space over 


( 2 ) 

(3) 


R, taking into account the boundary condition (1.8). 

the linear mapping Hp : Xx>p —> L 2 (Q) gives the reconstructed function, 
the linear mapping Vp : Xp i0 —> L 2 (tt) d gives a reconstructed gradient, 
which must be defined such that ||Vp • II^q^ is a norm on Xp^. 

Definition 2.6. (Gradient scheme for the obstacle problem). Let V be a gradient 
discretisation in the sense of Definition|2.5| The corresponding gradient scheme for 


( |2.2[ ) is given by 

Find u £ /Cp := {n £ Xp,o : npi> < g in fi} such that, Vu £ K. p, 

/ A(x)X7t>u(x) ■ Vp(u — v)(x) da; < / /(x)np(w — v)(x) da;. 

Jn Jo. 


(2.7) 
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The coercivity, consistency and limit-conformity of a gradient discretisation in 
the sense of Definition [23] are defined through the following constant and functions: 


(2.8) 

C v 

(2.9) 

V(v?, v) £ K x ICd, 

and 


(2.10) 


Vt/> £ 

H dW (n), WdW = 


V 


||n P u|| L 2 (n) 

max —n-> 

»6^D l0 \)0} || ^'DV\\L 2 (n) d 


St,{<p,v) = \\n. v v - p\\l 2 (Q) + \\Vw - V^|| L2(n)d , 


sup 

eA'c.oVfO} II V-D'C|| L 2( n) d 


/ (V'pv-'ip + IIxiV-divt'ip)) dx 

J n 


The only indicator that changes with respect to [Hi is St>■ For PDEs, 
the consistency requires to consider Sd on Hq(D) x Xd,o an( i to ensure that 
\ni ve x- Dm —> 0 as m —> oo for any ip £ Hq(D). Here, the domain of 

Sx> is adjusted to the set to which the solution of the variational inequality belongs 
(namely 1C), and to the set in which we can pick the test functions of the gradient 
scheme (namely ICd)- 

We note that this double reduction does not necessarily facilitate, with respect 
to the PDEs case, the proof of the consistency of the sequence of gradient discreti¬ 
sations. In practice, however, the proof developed for the PDE case also provides 
the consistency of gradient discretisations for variational inequalities. 


2.3. Error estimates. We present here our main error estimates for the gradient 


schemes approximations of Problems (2.1) and (2.2). 


2.3.1. Signorini problem. 


Theorem 2.7 (Error estimate for the Signorini problem). Under Hypothesis 2.1 


let u £ 1C be the solution to Problem \2. I\ j. If V is a gradient discretisation in the 
sense of Definition \2.S\ and ICd is non-empty, then there exists a unique solution 


u £ lC-D to the gradient scheme (2.3). Furthermore, this solution satisfies the 


following inequalities, for any vd £ ICd’- 
(2.11) ||Vx>u - Vu|| L 2( Q ),j < 


and 


^ Gd(u,vd) + + ^[Wx>(AVu) + (A + X)Sd(u, v v )\, 


(2.12) ||nx>u-u|| L 2 (n) < 


where 

(2.13) 


Cv\l ~^Gd{u, vd) + + -^{CdWd(HVu) + (Cd A + A)Sx> (u, vd)\, 


Gd(u,vd) = { r )n{XS7 u) ,TdVd — 7 u) and G jt, = max(0,G-p). 


Here, the quantities Cd, Sd(u, vd ) and Wd are defined by (2.4), (2.5) and (2.6). 
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2.3.2. Obstacle problem. An error estimate for gradient schemes for the obstacle 
problem is already given in [T, Theorem 1], For low-order methods with piecewise 
constant approximations, such as the HMM schemes, this theorem provides an 
0(y/h) rate of convergence for the function and the gradient (h is the mesh size). 
The following theorem improves [Tj Theorem 1] by introducing the free choice of 
interpolant up. This enables us, in Section |~2.4.2[ to establish much better rates of 
convergence - namely 0(h) for HMM, for example. 


Theorem 2.8 (Error estimate for the obstacle problem). Under Hypothesis 2.2 
let u £ K, be the solution to Problem 


_2.2j ). If V is a gradient discretisation in the 
sense of Definition \2.S[ and ifK -p is non-empty, then there exists a unique solution 
u £ K .-p to the gradient scheme (2.1). Moreover, if drv(A'Vu) £ L 2 (fl ) then this 
solution satisfies the following estimates, for any up £ IC'D- 


(2.14) ||Vdu-V tt|| L 2 (n) d < 


—Ed(u, up) + + t[Wb(AVu) + (A + A)iSp(u, up)], 
A A 


(2.15) lin^u - it|| L 2 ( n) < 

Cp \J-^Ex>(u,vx>) + + ^ [Cp Wp(AVu) + (CpA + A)iSp(u, up)], 

where 


(2.16) Ep(u,up) = / (div(AVii)+ /)(«, —Hpup) da; and E^, = max(0, Ed). 
an 

Here, the quantities Cp, 5p(u,up) and Wp are defined by (2.8|, (2.9) and (2.10|). 


Remark 2.9. We note that, if A is Lipschitz-continuous, the assumption div(AVu) £ 
L 2 (H) is reasonable given the l? 2 -regularity result on ii of |]. 


Compared with the previous general estimates, such as Falk’s Theorem [2S], our 
estimates seem to be simpler since it only involves one choice of interpolant up 
while Falk’s estimate depends on choices of Vh £ Kh and v £ K. Yet, our estimates 
provide the same final orders of convergence as Falk’s estimate. We also note that 
Estimates (2.11), (2.12), (2.14) and (2.15) are applicable to conforming and non- 
conforming methods, whereas the estimates in 361 seem to be applicable only to 
conforming methods. 


2.4. Orders of convergence. In what follows, we consider polytopal open sets 
and gradient discretisations based on polytopal meshes of D, as defined in mm- 
The following definition is a simplification (it does not include the vertices) of the 
definition in [191 . that will however be sufficient to our purpose. 

Definition 2.10 (Polytopal mesh). Let be a bounded polytopal open subset of 
M. d (d > 1). A polytopal mesh of H is given by T = (M,£,V), where: 

(1) A4 is a finite family of non empty connected polytopal open disjoint subsets 
of S2 (the cells) such that f l = UkgmK- F° r an y A' £ A4, \K\ > 0 is the 
measure of K and hx denotes the diameter of K. 
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(2) £ is a finite family of disjoint subsets of 12 (the edges of the mesh in 2D, 
the faces in 3D), such that any a £ £ is a non empty open subset of 
a hyperplane of and er C 12. We assume that for all K £ Ad there 
exists a subset £k of £ such that dK = U aG £ K a. We then denote by 
M<j = {K £ Ad : cr £ £k}- We then assume that, for all a £ £, Ad a has 
exactly one element and <r C <912, or Ad a has two elements and <r C 12. We 
let £; nt be the set of all interior faces, i.e. a £ £ such that er C 12, and £ ext 
the set of boundary faces, i.e. a £ £ such that a C <912. For a £ £, the 
{d— 1)-dimensional measure of a is |oj, the centre of gravity of a is x a , and 
the diameter of <r is h a . 

(3) V = (xk)kgM is a family of points of 12 indexed by Ad and such that, for 
all K £ Ad, xk € I\ (xk is sometimes called the “centre” of K). We then 
assume that all cells K £ Ad are strictly icif-star-shaped, meaning that if 
x is in the closure of K then the line segment [xk,x) is included in K. 

For a given K £ Ad, let n k,u be the unit vector normal to a outward to K and 
denote by dx,a the orthogonal distance between Xk and a £ £k■ The size of the 
discretisation is Hm = sup{hK ■ K £ Ad}. 

For most gradient discretisations for PDEs based on first order methods, includ¬ 
ing HMM methods and conforming or non-conforming finite elements methods, 
explicit estimates on St> and Wx> can be established m- The proofs of these es¬ 
timates are easily transferable to the above setting of gradient discretisations for 
variational inequalities, and give 

( 2 . 17 ) Wv(il>)<Ch M HUnny, ^£H\n) d , 


(2.18) mf Vr>e ^S v (ip,VT>) < Ch M \\<p\\H 3 (si), V<p £ i? 2 (12) n K. 

Hence, if A is Lipschitz-continuous and u £ H 2 (£L), then the terms involving Wz> 
and Sz> in ( |2.11| ) , ( |2.12[ ), ( |2.14[ ) and ( 2.15| ) are of order provided that vx> 

is chosen to optimise Sx>- Since the terms Gx> and Ed can be bounded above by 
(7117(11) — Tdvd\\l 2 (c)SI) and C\\u — n-pnx)|| l 2 (£ 1 ), respectively, we expect these to 

be of order Hm . Hence, the dominating terms in the error estimates are 


Gj and 

1 e£. In first approximation, these terms seem to behave as yfhM for a first order 
conforming or non-conforming method. This initial brute estimate can however 
usually be improved, as we will show in Theorems |2.11| and |2.13[ even for non- 
conforming methods based on piecewise constant reconstructed functions, and lead 
to the expected 0{Hm) global convergence rate. 


2.4.1. Signorini problem. For a first order conforming numerical method, the Pi 
finite element method for example, if u is in then the classical interpolant 

vd £ Xd,Fx constructed from the values of u at the vertices satisfies [7] 

IIHdcd — u||i2(n) < Ch 2 M \\D 2 u\\ L 2(£iyxd, 

||Vx>ud — Vu|| L 2 (n) d < Ch M \\D 2 u\\ L 2( n ' ) dxd, 

\\Tdvv ~ 7(w)IU 2 (9n) < Ch 2 M \\D 2 'y(u)\\ L 2( dn ' j dxd. 


If a is piecewise affine on the mesh then this interpolant vd lies in ICd and can there¬ 
fore be used in Theorem 2.7 We then have Sd(u,Vd) = 0(h M ) and Gx>(l 1,vd) = 
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0(h 2 M ), and (2.11 1 -(2.12) give an order one error estimate on the H 1 norm. If a is 


not linear, the definition of /Cp is usually relaxed, see Section [ 6 ] 

A number of low-order methods have piecewise constant approximations of the 
solution, e.g. finite volume methods or finite element methods with mass lumping. 
For those, there is no hope of finding an interpolant which gives an order h? M 
approximation of u in L 2 (Cl) norm. We can however prove, for such methods, 
that Gp behaves better than the expected order. In the following theorem, 

we denote by IF 2,00 (dCl) the functions v on dCl such that, for any face F of dCl, 
v\F £ W 2 ’°°{F). 

Theorem 2.11 (Signorini problem: order of convergence for non-conforming recon¬ 
structions). Under Hypothesis \2. 1\ let V be a gradient discretisation in the sense of 
Definition \2.3\ such that 1C p 0, and let T = (Ai,£,V) be a polytopal mesh of Cl. 
Let u and u be the respective solutions to Problems <H3 and (|2.3[). We assume that 


the barrier a is constant, that 7 n (AVh) £ L 2 (dCl) and that 7(h) £ W 2, °°(dCl). We 
also assume that there exists an interpolant vx> £ 1C p such that St>(u,vt>) < C\hM 
and that, for any a £ £ ext , ||TpPp - j(u)(x a )\\ L 2 ^ < C 2 h 2 |<r| 1 / 2 , where x„ £ a 
(here, the constants Ci do not depend on the edge or the mesh). Then, there exists 
C depending only on Cl, A, u, a, C \, C 2 and an upper bound of Gp such that 


(2.19) 


| Vp« ^ Vu||x,2(fj)<i + ||IIpU — li||L 2 (f2) < CHm + GWp(AVli). 


Remark 2.12. Since 7 „(AVw) £ L 2 (dCl) the reconstructed tr ace Tp can be taken 
with values in L 2 (dCl) (see the discussion at the end of Section 2.2.1). In particular, 
Tpu can be piecewise constant. The H l / 2 (dC£) norm in the definition (2.4) of Gp 
is then replaced with the L 2 (dCl) norm, and the duality product in (2.6) is replaced 
with a plain integral on dCl. 


For the two-point finite volume method, an error estimate similar to (2.19) is 
stated in [32], under the assumption that the solution is in H 2 (CL). It however seems 
to us that the proof in [32] uses an estimate which requires 7 (u) £ W 2 '°°(dCl), as 
in Theorem 12.Ill 

We notice that the assumption that a is constant is compatible with some models, 
such as an electrochemical reaction and friction mechanics [28]. This condition on 
a can be relaxed or, as detailed in Section [6j a can be approximated by a simpler 
barrier - the definition of which depends on the considered scheme. 


2.4.2. Obstacle problem. As explained in the introduction of this section, to esti¬ 
mate the order of convergence for the obstacle problem we only need to estimate 
Ex>(u, vx>). This can be readily done for Pi finite elements, for example. If g is 
linear or constant, letting rip = u at all vertices (as in M) shows that up is an 
element of 1C p, and that Sx>(u,v p) = 0(hj^ 1 ) and Et>(u,vt>) = 0(h 2 M ). Therefore, 
Theorem |2.8| provides an order one error estimate. 

The following theorem shows that, as for the Signorini problem, the error es¬ 
timate for the obstacle problem is of order one even for methods with piecewise 
constant reconstructed functions. 


Theorem 2.13 (Obstacle problem: order of convergence for non-conforming recon¬ 
structions). Under Hypothesis 2.2 let T> be a gradient discretisation in the sense of 
Definition\2.5\ such that 1C p 0, and let T = (Ai,£,V) be a polytopal mesh of Cl. 
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Let u and u be the respective solutions to Problems \2.2\ and ( |2.7[ ). We assume that 
g is constant, that div(AVu) £ L 2 (Tt), and that u £ We also assume 

that there exists an interpolant up £ /Cp such that St>(u,vt>) < C\hj ^ and, for any 
I\ £ A4, ||Hpup — u(xk)\ \l 2 (k) < C 2 h 2 K \K\ 1 / 2 (here, the various constants Ci do 
not depend on the cell or the mesh). Then, there exists C depending only on f2, A, 
A, Ci, C 2 and an upper bound o/Cp such that 

(2.20) ||Vpu — Vu|| iW + ||IIpu — w||l2(q) < ChM+CWv{AS7u). 


As for the Signorini problem, under regularity assumptions the condition that 
g is constant can be relaxed, see Remark |5.1| However, if g is not constant it 
might be more practical to use an approximation of this barrier in the numerical 
discretisation - see Section [ii] 


Remark 2.14 (Convergence without regularity assumptions). The various regularity 
assumptions made on the data or the solutions in all previous theorems are used 
to state optimal orders of convergence. The convergence of gradient schemes can 
however be established under the sole regularity assumptions stated in Hypothesis 

see [I]. 


2.1 and 2.2 


3. Examples of gradient schemes 


Many numerical schemes can be seen as gradient schemes. We show here that 
two families of methods, namely conforming Galerkin methods and HMM methods, 
can be recast as gradient schemes when applied to variational inequalities. 

3.1. Galerkin methods. Any Galerkin method, including conforming finite ele¬ 
ments and spectral methods of any order, fits into the gradient scheme framework. 

For the Signorini problem, if V is a finite dimensional subspace of {z> £ H 1 (fl) : 
j(v) = 0 on Ti}, we let Any, = V, Hp = Id, Vp = V and Tp = 7. Then 
prp,n,np,Tp,Vp) is a gradient discretisation and the corresponding gradient 


scheme corresponds to the Galerkin approximation of (2.2). Then, Cp defined by 


(2.4) is bounded above by the maximum between the Poincare constant and the 
norm of the trace operator 7 : H 1 (LI) -A iJ 1 / 2 (9S2). 

Conforming Galerkin methods for the obstacle problem consist in taking Apy = 
V a finite dimensional subspace of Hg(fl), Hp = Id and Vp = V. For this gradient 
discretisation, Cp (defined by (2.81) is bounded above by the Poincare constant in 

H^n). 

For both the Signorini and the obstacle problems, Wp is identically zero and the 
error estimate is solely dictated by the interpolation error 5p(u, v) and by Gp(S, v) 
or Ex>(u,v), as expected. For Pi finite elements, if the barrier (a or g) is piecewise 
affine on the mesh then the classical interpolant v constructed from the nodal values 


of u belongs to 1C. As we saw, using this interpolant in Estimates (2.14) and (2.15) 
leads to the expected order 1 convergence. 

If the barrier is more complex, then it is usual to consider some piecewise affine 
approximation of it in the definition of the scheme, and the gradient scheme is then 
modified to take into account this approximate barrier (see Section[6]); the previous 
natural interpolant of u then belongs to the (approximate) discrete set /Cp used to 
define the scheme, and gives a proper estimate of Sx>{u, v) and Gti{u, v) or £p(zt, v). 
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3.2. Hybrid mimetic mixed methods (HMM). HMM is a family of methods, 
introduced in [16] , that contains the hybrid finite volume methods |23| , the (mixed- 
hybrid) mimetic finite differences methods [5] and the mixed finite volume methods 
M- These methods were developed for anisotropic heterogeneous diffusion equa¬ 
tions on generic grids, as often encountered in engineering problems related to flows 
in porous media. So far, HMM schemes have only be considered for PDEs. The 
generic setting of gradient schemes enable us to develop and analyse HMM methods 
for variational inequalities. 

Remark 3.1. The (conforming) mimetic methods studied in J2] for obstacle problems 
are different from the mixed-hybrid mimetic method of the HMM family; they are 
based on nodal unknowns rather than cell and face unknowns. These conforming 
mimetic methods however also fit into the gradient scheme framework [TUI . 

3.2.1. Signorini problem. Let T be a polytopal mesh of Q. The discrete space to 
consider is 

Xv.Ti = = ((vk)kgM, ('Va)aGs) ■ Vk G R, V a G R, 

v a = 0 for all a G £ ext such that cr C T^ 

(we assume here that the mesh is compatible with (Pi)j=i 2 , 3 , in the sense that for 
any i = 1,2,3, each boundary edge is either fully included in T^ or disjoint from 
this set). The operators Hp and Vp are defined as follows: 

\/v G Xp^! , VA' G M. : n v v = vk on AT, 


Vu G Xp >ri , VAT G M : Vpw = V K v + ^-{A K R K {v)) a n K ,a on co({cr, x K }) 


l-K,, 

where co (S) is the convex hull of the set S and 

• v kv = J2 a& e K W\(v a - v K )n K>cr , 

• Rk(v) = (v„ - V K ~ Vjs'O • ( X a - X K ))<t££ k G 

• Ak is an isomorphism of the vector space Im(A^-). 

It is natural to consider a piecewise constant trace reconstruction: 


I!ard (£k) 


(3.2) 


Vcr G fext : Tpi; = v„ on a. 


1/2 

This reconstructed trace operator does not take values in H r (<9fl), and the cor¬ 
responding gradient discretisation T> is therefore not admissible in the sense of 
Definition |2.3| However, for sufficiently regular u, we can consider reconstructed 


2.12 


traces in L 2 (dfl) only, see Remark 

It is proved in m that, if Ti = dfi, then (Xpy^np, Vp) is a gradient dis¬ 
cretisation for homogeneous Dirichlet boundary conditions, whose gradient scheme 
gives the HMM method when applied to elliptic PDEs. 

With /Cp := {u G Xpy, : v a < a on a , for all a G £ ext such that cr C Ta}, the 
gradient scheme (2.3) corresponding to this gradient discretisation can be recast as 


(3.3) 


Find u G K. p such that, for all v G K. p , 

|A \Ak\7ku ■ Vx [u — v) + 'y ' Rk(u — v) t M k R k (u) 

k^m 


K&M 


< ^2 ( u k~v K ) / /( x)dx. 

KGM J k 
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where A k is the value of A in the cell K (it is usual - although not mandatory 
- to assume that A is piecewise constant on the mesh) and, for K £ Ad, is a 
symmetric positive definite matrix of size Card (£k)- This matrix is associated to 
Ak, see [IS] for the details. 

Under classical assumptions on the mesh and on the matrices Mk, it is shown 
in JT8] that the constant Cp can be bounded above by quantities not depending 
on the mesh size, and that Wp(i/>) is of order O(hj^). If d, < 3 and u £ H 2 (£l), 
we can construct the interpolant up = ((vk)k<eMi ( v cr)aeM) with vk = u(xk) and 
v a = u{x,y ) and, by m Propositions 13.14 and 13.15], we have 


Sx>(u,Vt>) — ||HpUp — |z, 2 (f2) + ll^pUp ~ Vu|| L 2( n )d < Ch M \\u\\ H 2^y 

Under the assumption that u £ 1U 2,00 (U), this estimate is also proved in [TO] (with 
INIir 2 (n) replaced with 11|(r*))- If u £ K. and a is constant, this interpolant 


vt> belongs to /Cp. Moreover, given the definition (3.21 of Tp, we have (Tpup)^ — 
7 (u)(x a ) = v a — u(x cr ) = 0. Hence, using this up in Theorem 2.11 we see that the 
HMM scheme for the Signorini problem enjoys an order 1 rate of convergence. A 
non-constant barrier a can be approximated by a piecewise constant barrier on £ ext 


and Theorem 6.2 shows that, with this approximation, we still have an order 1 rate 
of convergence. Note that since this convergence also involves the gradients, a first 
order convergence is optimal for a low-order method such as the HMM method. 


3.2.2. Obstacle problem. We still take T a potytopal mesh of f2, and we consider 
V = (Xp, 0 ,n p , Vp) given by Xt>,o = Xx>,an (see dO)) and Hp and Vp as in 


Section 3.2.1 


Using V in (2.7), we obtain the HMM method for the obstacle 
vk < g on K, for all K £ A4}, this method 


problem. Setting K ,p := {r> £ Xp 0 
reads 


(3.4) 


Find u £ /Cp such that, for all v £ /Cp , 

\K\A k Vku ■ Vk(« — v) + ^ Rk(u — v) t M k R k (u) 


kgm 


kgm 


< ^2 ( u k~vk) / f(x) dx 

K£M 


Using the same interpolant as in Section 3.2.1 we see that Theorem 2.13 (for a 


constant barrier g) and Theorem 6.4 (for a piecewise constant approximation of g) 


provide an order 1 convergence rate. 

To our knowledge, the HMM method has never been considered before for the 
Signorini or the obstacle problem. These methods are however fully relevant for 
these models, since HMM schemes have been developed to deal with anisotropic 
heterogeneous diffusion processes in porous media (processes involved in seepage 
problems for example). The error estimates obtained here on the HMM method, 
through the development of a gradient scheme framework for variational inequali¬ 
ties, are therefore new and demonstrate that this framework has applications be¬ 
yond merely giving back known results for methods such as Pi finite elements. 


4. Numerical results 


We present numerical experiments to highlight the efficiency of the HMM meth¬ 
ods for variational inequalities, and to verify our theoretical results. To solve the 
inequalities (3.4) and (3.3), we use the monotony algorithm of [50s. We introduce 
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the fluxes (Fx,o(u))x^M,crG£ K defined for u £ FIp {Hp = Xp t r 1 for the Signorini 
problem, Hp = Xp$ for the obstacle problem) through the formula [HJ Eq. (2.25)] 


VK £ M, Vv £ Hp : 

y, |cr| F k ,o{u){v k - V„) = \K\A k X7 k u ■ S7 K V + Rx(u)M K R K (v). 

cr&Sl< 


Note that the flux F Ka thus constructed is an approximation of — AV u . n K a . 

The variational inequalities (3.4) and (3.3) can then be re-written in terms of the 
balance and conservativity of the fluxes in the following ways: 


• Signorini problem: 


y W\F k A u ) = m(K)f K , 

<j££ k 

FrA*) + f lA u ) = 

U(j — 0, 

FkA u ) = 

FK,a{F)(U(j Oj(j^ — 0 , 

-FkA u ) < 0, 

U& — O'er j 


\/KeM , 

Vcr G £ int with M a = {K , L}, 

V(T £ fext such that a C Id, 

VAT £ M ,\/a £ £k such that cr C r 2 , 
Vif £ M ,Va £ £k such that cr C T 3 , 
Vif £ M ,\/a £ £r such that cr C T 3 , 
Vcr £ fext such that cr C f 3 . 


Here, a a is a constant approximation of a on cr. 
• Obstacle problem: 



cr F K ,oi u ) + ( 9K - U K ) = o, 

\/K £ M, 

/ 

F k ,o(u) + FlAu) = o, 

Vcr £ 

y M f k A u ) < A k )/k, 

VK £ M, 

ctGSk 


UK < gK, 

VK £ M, 

U(j — 0, 

Vcr £ £ e xt • 


with Mo = {K,L}> 


Here, gx is an approximation of g on K , see Section [6] 

The iteration monotony algorithm detailed in |30| can be naturally applied to 
these formulations. We notice that by |371 Section 4.2] the HMM method is mono¬ 
tone when applied to isotropic diffusion on meshes made of acute triangles; in that 
case, the convergence of the monotony algorithm can be established as in [30] for 
finite element and two-point finite volume methods. In our numerical tests, we 
noticed that, even when applied on meshes for which the monotony of HMM is 
unknown or fails, the monotony algorithm actually still converges. 

All tests given below are performed by MATLAB code using the PDE toolbox. 
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Figure 4.1. The HHM solution for Test 1 with h = 0.05. The 
line represents where the boundary condition on T 3 changes from 
Neumann to Dirichlet. 


Test 1. We investigate the Signorini problem from [40] : 


—A-u = 27 rsin( 27 ra:) 

u = 0 

vs -n = o 


in Cl = (0, l) 2 , 
on Ti, 
on r 2 , 


«>0 | 

Vu • n > 0 > on T 3 , 
uS7u • n = 0 J 


with Ti = [0,1] x {1}, T 2 = ({0} x [0,1]) U ({1} x [0,1]) and T 3 = [0,1] x {0}. Here, 
the domain is meshed with triangles produced by INITMESH. 

Figure [~i~T| presents the graph of the approximate solution obtained on a mesh 
of size 0.05. The solution compares very well with linear finite elements solution 
from |42| : the graph seems to perfectly capture the point where the condition on 
r 3 changes from Dirichlet to Neumann. Table [l] shows the number of iterations 
(NITER) of the monotony algorithm, required to obtain the HMM solution for var¬ 
ious mesh sizes. Herbin in [321 proves that this number of iterations is theoretically 
bounded by the number of edges included in T 3 . We observe that NITER is much 
less than this worst-case bound. We also notice the robustness of this monotony 
algorithm: reducing the mesh size does not significantly affect the number of iter¬ 
ations. 










16 


YAHYA ALNASHRI AND JEROME DRONIOU 


Table 1. Number of iterations of the monotony algorithm in Test 1 


Mesh size h 

0.0625 

0.0500 

0.0250 

0.0156 

#{cr c r 3 } 

16 

20 

40 

46 

NITER 

4 

5 

5 

6 




Figure 4.2. Test 2: geometry and diffusion (left), typical grid (right). 


Test 2. Most studies on variational inequalities, for instance [39]|40], consider test 
cases for which the exact solutions are not available; rates of convergence are then 
assessed by comparing coarse approximate solutions with a reference approximate 
solution, computed on a fine mesh. To more rigorously assess the rates of con¬ 
vergence, we develop here a new heterogeneous test case for Signorini boundary 
conditions, which has an analytical solution with non-trivial one-sided conditions 
on T3 (the analytical solution switches from homogeneous Dirichlet to homogeneous 
Neumann at the mid-point of this boundary). 

In this test case, we consider here PI p| with the geometry of the domain 
presented in Figure [42} left. 

The exact solution is 


f P(y)h(x) for y G (0, \) and x G (0,1), 
\ xg(x)G(y) for y G (|, 1) and z G (0,1). 


To ensures that (1.21 holds on the lower part of Ti, and that the normal derivatives 
AVtt-n at the interface y = \ match, we take P such that P(0) = P{\) = P'{\) = 0. 
Assuming that P < 0 on (0, |), the conditions d n u = —d x u = 0 and u < 0 on the 
lower half of T3, as well as the homogeneous condition d n ti = 0 on the lower half 
of r2, will be satisfied if h is such that h'{ 0) = h'( 1) = 0 and h{ 0) = 1. We choose 



and h(x) = cos(nx). 


Given our choice of P , ensuring that div(AVu) G L 2 (I2) (i.e. that u and the normal 
derivatives match at y = ^) c an be done by taking G such that G(^) =G'{\) =0. 
The boundary condition (1.2) on the upper part of Fi is satisfied if G(l) = 0. The 
boundary conditions d n ii = 0 on the upper part of r 2 , and u = 0 and d n u < 0 
on the upper part of r 3 , are enforced by taking g such that g( 1) + </(l) = 0 and 
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ff(0) > 0 (with G > 0 on (|, 1)). We select here 
' 7r \ l+cos(7rx) 


g{x) = cos 2 = 


and G{y) = (1 - y) ( y - - 


This u is then the analytical solution to the following problem 



—div(AVw) = / 

in f2, 


u = 0 

on Ti, 


AVu ■ n = 0 

on r 2 , 


u < 0 ) 
AVm • n < 0 > 
uAVu • n = 0 J 

on r 3 , 

with, for y < 



f{x,y) = 100 

n 2 y cos(nx) ^y — — 2ycos(nx) 

and, for y > 




f{x,y) = — cos(Trx)(y — 1) (y- -2a; cos 2 {3y - 2) 


+ 7rsin(7rx)(y- 1) ( y - - 



Figure 4.3. The HHM solution for Test 2 on an hexagonal mesh 
with h = 0.07. 

We test the scheme on a sequence of meshes (mosly) made of hexagonal cells. 
The third mesh in this sequence (with mesh size h = 0.07) is represented in Figure 
|4.2[ right. 

The relative errors on u and V«, and the corresponding orders of convergence 
(computed from one mesh to the next one), are presented in Table [2j For the 
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HMM method, Theorem |2.11| predicts a convergence of order 1 on the gradient 
if the solution belongs to H 2 and AVu £ H 1 ; the observed numerical rates are 
not slightly above these values, probably due to the fact that A is not Lipschitz- 
continuous here (and thus the regularity AVi2 £ H 1 is not satisfied). As in the 
previous test, the number of iterations of the monotony algorithm remains well 
below the theoretical bound. 

Table 2. Error estimate and number of iterations for Test 2 


mesh size h 

0.24 

0.13 

0.07 

0.03 

ll“-n-DU|| t 2 (!2) 

Order of convergence 

0.6858 

0.2531 

1.60 

0.1355 

0.92 

0.0758 

0.84 

!|V'U-Vx>'!i|| L 2 (SJ) 

lt- v “IK2 (n) 

Order of convergence 

0.4360 

0.2038 

1.22 

0.1041 

0.99 

0.0542 

0.95 

#{<t c r 3 } 

20 

40 

80 

160 

NITER 

5 

5 

6 

7 


The solution to the HMM scheme for the third mesh in the family used in Table 
[2] is plotted in Figure |4~3) On T 3 , the saturated constraint for the exact solution 
changes from Neumann VS • n = 0 to Dirichlet St = 0 at y = 0.5. It is clear 


on Figure 4.3 that the HMM scheme captures well this change of constraint. The 
slight bump visible at y = 1/2 is most probably due to the fact that the mesh is not 
aligned with the heterogeneity of A (hence, in some cells the diffusion tensor takes 
two -very distinct- values, and its approximation by one constant value smears the 
solution). 

When meshes are aligned with data heterogeneities, such bumps do not appear. 
This is illustrated in Figure |4~4} in which we represent the solution obtained when 
using a “Kershaw” mesh as in the FVCA5 benchmark [31]. This mesh has size 
h = 0.16, 34 edges on r 3 , and the monotony algorithm converges in 7 iterations. 
The relative L 2 error on u and VS are respectively 0.017 and 0.019. As expected 
on these kinds of extremely distorted meshes, the solution has internal oscillations, 
but is otherwise qualitatively good. In particular, despite distorted cells near the 
boundary F 3 , the transition between the Dirichlet and the Neumann boundary 
conditions is well captured. 


Test 3. In this case test, we apply the HMM method to the homogeneous obstacle 
problem: 

(Au + C)(ip — u) = 0, 

-A« > C, 
u> ip, 
u = 0, 

The constant C is negative, and the obstacle function ip(x , y) = — min(x, 1— x, y , 1— 
y) satisfies ip = 0 on 9D. 


in fl = (0, l) 2 , 
in Q, 
in n, 
on dfl. 


Figure 4.5 shows the graph of the HMM solution to the above obstacle problem 
with respect to h = 0.05. Table [3] illustrates the performance of the method and the 
algorithm. Here again, the number of iterations required to obtain the solution is 
far less than the number of cells, which is a theoretical bound in the case of obstacle 
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Figure 4.4. The HHM solution (left) for Test 2 on a Kershaw 
mesh (right). 


problem [3Uj. Our results compare well with the results obtained by semi-iterative 
Newton-type methods in [5], which indicate that decreasing \C\ contributes to the 
difficulty of the problem (leading to an increase number of iterations). We note 
that, for a mesh of nearly 14,000 cells, we only require 29 iterations if \C\ = 5 and 
14 iterations if \C\ = 20. On a mesh of 10,000 cells, the semi-iterative Newton-type 
method of [9) requires 32 iterations if \C\ = 5 and 9 iterations if \C\ = 20. 



Figure 4.5. The HHM solution for Test 3 with h = 0.05, C = —20. 
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Table 3. The number of iterations for various C 


mesh size 

h 

0.016 

0.025 

0.050 

0.062 



14006 

5422 

1342 

872 

NITER (C = 

= -5) 

29 

20 

12 

11 

NITER (C = 

-10) 

23 

18 

10 

9 

NITER (C = 

-15) 

19 

13 

9 

8 

NITER (C = 

-20) 

14 

12 

7 

8 


5. Proofs of the theorems 


Proof of Theorem |2.7[ Since /Cp is a non-empty convex closed set, the existence 
and uniqueness of a solution to Problem (2.3) follows from Stampacchia’s theorem. 
We note that AVzi £ H^ v (fl) since 


(5.1) 


— div(AVu) = / in the sense of distributions 


(use v = u + xp in (2.1), with ip £ C'^°(f2)). Taking xp = AVu in the definition (2.6) 
of Wt> , for any z £ Ic,rt we get 

(5.2) / Vpz-AVudx-l- / IIp,s • div(AVw) dx < 

in in 

II Vp^|| i 2(Q)dWp(AVu) + (7 n (AVu), Tpz). 

Let us focus on the last term on the right-hand side of the above inequality. For 
any up € /Cp, we have 


(5-3) (7„(AVu), Tp(up - u)) = 


(7n(AVu),Tpup - 7u) + (7„(AVu),7ii - Tpit). 

The definition of the space shows that there exists w £ H 1 ( fl) such that 

7ru = Tpit (this implies that w £ K.). According to the definition of the duality 
product dO, we have 


(7n (.AV u), 'yu — Tpit) = / AVu, • V(zt — w) da; + / div(AVu)(u — u>) dx. 
in in 

Since u satisfies (5.1), it follows that 

(5.4) ( 7n (AVu),7u-Tpu) = / AVu • V(zt — ru) dx — / f(u — w) dx, 

in in 


which leads to, taking re as a test function in the weak formulation (2.1), 

(5.5) (7n(AVu), yu - Tpu) < 0. 


From this inequality, recalling (5.1) and (5.3) and setting z = up — u £ Xp in 
(5.2), we deduce that 

(5.6) / Vp(up — u) • AVu dx + / /IIp(u — up) dx < 

in in 

II Vp(up — u)|| L 2( n )dWp(AVu) + Gp(u, up). 
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We use the fact that u is the solution to Problem (2.31 to bound from below the 
term involving / and we get 

/ AVp(i?p — u) ■ (Vzi — Vptt) da? < 

J n 

|| Vp(up — u)|| i 2(f 2 )dW : D(AVu) + Gp(u, vt)) + 

(note that we also used the bound Gp < G)(,). Adding and subtracting Vpup in 
Vu — \7tiu, and using Cauchy-Scliwarz’s inequality, this leads to 

AIIVpup - Vpu||2 2(n)d < 

|| Vprp — Vpu||i2(Q)d (^Wt>{AS7u) + ASp(u, up)^ + Gp(u, up) + . 
Applying Young’s inequality gives 

_ [2 1 A 2 

(5.7) || Vpup — Vpu|| L 2 (Q)d < , / — Gp(u, up)+ + ^ [Wp(AVu) + A<Sp(«, up)] . 


Estimate (2.11 1 follows from || Vpit — Vu\\ L 2 (Q)d < ||Vpu—VpUp||p2(Q)d + ||VpUp — 
Vu|| i 2 (p)d and y/a + b < y/a + Vb. Using (2.4) and (|5.7|, we obtain 


2 1 _ 2 

HIIpup — IIpu||i2(Q) < Gp. / -Gp(S, Up)+ + [Wp(AVu) + A 5 p(zt, up)] . 

By writing ||IIpzt - u\\ L ^(U) < ||IIpu - npUp|| L 2 (n ) + ||IIpi;p - u||l2 ( q), the above 
inequality shows that 

l|Bpw - u||l2(q) < 


Gpy —Gp(u,up)+ + -^2 ^Wp(AVit) + ASp(u,up)^ + 5p(u,up)■ 

Applying y/a + b < y/a + y/b again, Estimate (2.12) is obtained and the proof is 
complete. 

Proof of Theorem 12.81 

The proof of this theorem is very similar as the proof of [1| Theorem 1], We 
however give some details for the sake of completeness. 

As for the Signorini problem, the existence and uniqueness of the solution to 
Problem \2.7\ follows from Stampacchia’s theorem. Let us now establish the error 
estimates. Under the assumption that div(AVii) £ L 2 (U), we note that AV« £ 
-Hdiv(^)- For any w £ Xp j0 , using ip = A Vu in the definition (2.10) of Wp 
therefore implies 

(5.8) / Vpw • AVudx + / np» • div(AVu) dx < ||Vpw|| L 2 (n) d Wp(AVu). 

J n J n 

For any up £ /Cp, one has 

(5.9) / Hp(it — up)div(AVu) da? = j (Ilpit — g)(div(AVu) +/) dx 

Jn Jo, 


+ f {g - npUp)(div(AVw) + /) da: - f /Hp(u - up) dx. 
J O JQ, 
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It is well known that the solution to the weak formulation (2.21 satisfies (1.6) in 


the sense of distributions (use test functions v = u — <p in (2.2), with ip £ C£°(fl) 


non-negative). Hence, under our regularity assumptions, (1.6) holds a.e. and, since 

u € K-ti, we obtain / (Hp u — g)(div(AVu) + /) da; < 0. Hence, 
an 


/ np(u-up)div(AVu) da; 

In 

< [ {g - npup)(div(AVu) + /) dx - f /Hp(u - up) da; 

Jn an 

= f (g — tt)(div(AVu) + /) da; + f [u — Hpi>p)(div(AVu) + /) da; 
an an 

- / /np(u - vd) dx. 
an 


Our regularity assumptions ensure that u satisfies (1.5). Therefore, by definition of 

Ed, 


(5.10) / ELd{u — up)div(AVu) da; < Ed{u, vd) — / /%(« — vd) dx. 

J Q J Q 


From this inequality and setting w = vd — u £ Xp o i n (5.8), we obtain 


/ Vp(up — u) ■ AVatda; + / /Hd{u — vd) da; < 
an an 


||Vp(up — w)|| L 2(q)<jWp(AVu) + Ed(u, vd)- 


The rest of proof can be handled in much the same way as proof of Theorem |2.7| 


taking over the reasoning from (5.6). See also |T) for more details. 


Proof of Theorem 


2.11 


We follow the same technique used in [22. According 
to Remark [2.12| and due to the regularity on the solution, since Tpup = 7 (u) = 0 


on Ti and 7 n (AVu) = 0 on f2, we can write 


Gd(u, vd) 


(Tdvd - 7(w))7n(AVu) da;. 


We notice first that the assumptions on u ensure that it is a solution of the 
Signorini problem in the strong sense. Applying Theorem |2.7| we have 

(5.11) || Vpu — Vu|| i 2(Q)d + HHr, it — u||z,2(n) 

< C (Wd(EVu) + Sd(u,vd) + \/Gd(u, J'x>) + ^ 


with C depending only on A, A and an upper bound of Cp. Since Sp(it, up) < CiIt-m 
by assumption, it remains to estimate the last term Gp(zt,up) + . We start by 
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writing 


Gt>(u,vt>) = / - a)y n (AVw) dx + / (a - 7(w))y n (AVu) Ax 

J r 3 J r 3 

= / (Td^x) - a)y n (AVu) dx 
Jr 3 

= 5Z /(TdVj, - a)7 n (AV«) dx 

aG£, aCr 3 “' <T 

= : G a {u,V V ), 

cr&£ , crcr 3 


where the term involving (a — y(tt))y n (AVM) has been eliminated by using (1.4). 
We then split the study on each er depending on the cases: either (i) yu < a a.e. on 
er, or (ii) meas({y £ er : "fu(y) —a = 0}) > 0 (where meas is the {d— l)-dimensional 
measure). 


In Case (i), we have y n (AVw) = 0 in er since yu satisfies (1.4). Hence, G a (u, vx>) = 


0. 

In Case (ii), let us denote V CT the tangential gradient to er, and let us recall that, 
as a consequence of Stampacchia’s lemma, if w £ W 1,1 (a) then V CT w = 0 a.e. on 
{y £ a : w{y) = 0} (here, “a.e.” is for the measure meas). Hence, with w = yu — a 
we obtain at least one yo £ cr such that (y u — a)(yo) = 0 and \7 (J ('yu — a)(yo) = 0. 
Let F be the face of <9f2 that contains a. Using a Taylor’s expansion along a path 
on F between yo and x a , we deduce 




(5.12) |(yu - a)(x<j)\ < L F h 2 

where L F depends on F and on the Lipschitz constant of the tangential derivative 
Van(7M — a) on F (the Lipschitz-continuity of Van (yzi — a) on this face follows from 
our assumption that y u — a £ W 2,00 (9H)). Recalling that \\ r ^T>VT> — 'yu(x a )\\L 2 (a) < 
C 2 h 2 lo) 1 / 2 , we infer that ||T -pv-p — a\\ L i^ < Ch 2 M lerl 1 / 2 , and therefore that 
\G a (u,v-D)\ < C'/i^,||y„(AVit)|| i 2 ((T) |(7| 1 / 2 . Here, C depends only on !!, « and 
a. 

Gathering the upper bounds on each G a and using the Cauchy-Schwarz inequal¬ 
ity gives Gt>{u,v-d ) < | |y n (AVi2) | |z,2(ar2) ? with C depending only on H, A, u 

and a. This implies Gx>(u, vx>) + < Ch 2 M , and the proof is complete by plugging 
this estimate into (5.11). 


Proof of Theorem |2.13[ The proof is very similar to the proof of Theorem 
2.11 As in this previous proof, we only have to estimate Ex>(u,vx>), which can be 
re-written as 

E-p(u, v-d) = / (div(AVw) + f)(u — g)&x + / (div(AVu) + f)(g - Ux,vt>) da; 

Jn Jn 

= [ (div(AVrt) + f){g — UtiVv) da; 

Jn 

= 5Z / (div(k'Vu) + f){g-Ux>VT>)dx =: 5^ E k (u,vt>), 

TS r- K A ” K Ty s— \ A 


KGM 


kgm 


where the term involving (div(AVu) + f)(u — g ) has been eliminated using (1.5). 


Each term E F (u, v-p) is then estimated by considering two cases, namely: (i) either 
u < g on K , in which case Ek(u,Vt>) = 0 since / + div(Au) = 0 on K, or (ii) 
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|{y £ K : u(y) = g}| > 0, in which case Ek{u,vd) is estimated by using a Taylor 
expansion and the assumption ||IIpup — u(xk )| \l 2 (k) < C 2 h 2 K \I\\ 1 l 2 . 

Remark 5.1. If (xk)kgm are the centres of gravity of the cells, the assumption 
that g is constant can be relaxed under additional regularity hypotheses. Namely, 
if F := div(AVu) + / € H l (K) and g £ H 2 (K) for any cell K, then by letting Fk 
and gK be the mean values on K of F and g , respectively, we can write 

E k (u, up) = / F{g - Up up) dx 
J I< 

= [ F{g{x K ) ~ IIpUp)dx + [ (F - F K )(g - g K )dx 
Jk Jk 

+ [ f (9k - g{xK)) dx =: T 1>k + T 2 ,k + T 3! k. 

Jk 

The term T\ k is estimated as in the proof (since u{xk) — g{xx) can be estimated 
using a Taylor expansion about y o). We estimate T% k by using the Cauchy-Schwarz 
inequality and classical estimates between an H 1 function and its average: 

\ t 2 ,k\ < H^ 1 - F K \\ L 2 ( K )\\g - gxWmK) < Ch 2 K \\F\\ H i( K )\\g\\ H i^ K y 

As for T 3 K , we use the fact that g £ H 2 (VL) and that Xk is the center of gravity 
of K to write | g{xK) ~ 9 k\ < Ch 2 K \\g\\ H 2 ^ K y Combining all these estimates to 
bound Ek(u,vd) and using Cauchy-Schwarz inequalities leads to an upper bound 
inO(^)fo;^,u D ). 


6. The case of approximate barriers 

For general barrier functions (a in the Signorini problem, g in the obstacle prob¬ 
lem), it might be challenging to find a proper approximation of u inside ICx>- Con¬ 
sider for example the Pi finite element method; the approximation is usually con¬ 
structed using the values of u at the nodes of the mesh, which only ensures that 
this approximation is bounded above by the barrier at these nodes, not necessarily 
everywhere else in the domain. It is therefore usual to relax the upper bound im¬ 
posed on the solution to the numerical scheme, and to consider only approximate 
barriers in the schemes. 


6.1. The Signorini problem. We consider the gradient scheme (2.3) with, instead 
of IC-D, the following convex set: 


( 6 . 1 ) 


A-p = {u £ Xx>,r 1 : T dv < o-d on Ts}, 


where ax> £ L 2 (dfl) is an approximation of a. The following theorems state error 
estimates for this modified gradient scheme. 


Theorem 6.1 (Error estimates for the S igno rini problem with approximate bar- 

if K.t> is not empty then there exists 


2.7 


rier). Under the assumptions of Theorem 
a unique solution u to the gradient scheme (|2.3|) in which /C-p has been replaced 
withK-D- Moreover , if a — op £ H Vi {dCl), Estimates (2.11) and (2.12) hold for 
any up £ K. p, provided that Gp(u, up) is replaced with 


Gp(u,up) = Gp(u,up) + (7 n (AVu),a - ap) 
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Proof. The proof is identical to that of Theorem |2.7[ provided we can control the 


left-hand side of (5.3). For any up £ /Cp, we write 


(7n(AVu),Tp(up - u)) = (7„(AVu), Tpup - ju) 

+ ( 7 n (A Vu), 7 u - (Tpu + a- a®)) + ( 7 „(AVm), a - ap). 

We note that the first term in the right-hand side is Gp(it, up); hence, the first 
and last terms in this right-hand side sums up to Gp(u, up) and we only need 
to prove that the second term is non-positive. We take w £ such that 

7 (w) = Tpit + a — ap, and we notice that w belongs to /C since Tpu and a — ap 
both belong to iJj/ (9f2), and since Tpu + a — ap < a on T 3 . Similarly to (5.4), 
the definition of w shows that 

( 7 n(AVu), 7 'it - (Tpit + a - ap)) = / AVu • V(u — w) dx — / f(u — w)dx. 

Jn Jn 

Hence, by B ( 7 n(AVu), 7 u— (Tpu + a —ap)) < 0 and, as required, the left-hand 
side of (5.3) is bounded above by Gp(w,up). □ 


In the case of the Pi finite element method, ap is the Pi approximation of a on 
912 constructed from its values at the nodes. The interpolant up of u mentioned in 
Section ra is bounded above by a = ap at the nodes, and there fore everywhere 


by ap; vt> therefore belongs to /Cp and can be used in Theorem 


6.1 


Under H 2 


regularity of a, we have ||ap — a||p 2 (as 7 ) = Q{h 2 M ). Hence, if 7 n (AVu) £ L 2 (9f2), 


we see that Gp(u, up) = 0(h 2 M ) and Theorem 6.1 thus gives an order one estimate. 


For several low-order methods (e.g. HMM, see [18]), the interpolant up of the 
exact solution u is constructed such that Tpup = "f(u)(x a ) on each er £ £ e xt- The 
natural approximate barrier is then piecewise constant, and an order one error 
estimate can be obtained as shown in the next result. 

Theorem 6.2 (Signorini problem: order of convergence for piecewise constant 


reconstructions). Under Hypothesis 2.1 let T> be a gradient discretisation, in the 
sense of Definition 2.3 and let T = (M,£,V) be a polytopal mesh of 12. Let 
ap £ L 2 (9f2) be such that ap —a = 0 onT\ and, for any a £ £ ex t such that a C T 3 , 
ap = a(X( 7 ) on a, where x„ £ a. Let ii be the solution to (2.1), and let us assume 
that 7 n (AVu) £ L 2 (dfl) and that 7 (u) — a £ W 2,oc (dfl). We also assume that 
there exists an interpolant up £ Xpy, such that St>(u,Vt>) < C\hj^ with C\ not 
depending on V or T, and, for any a C r 3 , Tpt>p = 7 (u)(x a ) on a. 

Then /Cp 7 ^ 0 and, if u is the solution to ( |2.3| ) with /Cp replaced with /Cp, it 
holds 


(6.2) ||Vpu — Va|| i 2(Q)d + ||Hpu — u\\l^(q.) < Ch.M + GWp(AVu) 

where C depends only on 12, A, u, a, Gi and an upper bound of Cxi- 

Proof. Clearly, up £ /Cp. The conclusion follows from T heore m 
that Gp(u,up) + = 0{h 2 M ). Note that, following Remark 
consider a piecewise constant reconstructed trace Tpu. 


2.12 


6.1 


if we prove 
it makes sense to 
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With our choices of ap and Hpxp, and using the fact that u is the solution to 


the Signorini problem in the strong sense (and satisfies therefore (1.4)), we have 


Gp(u,xp) = / (Tpxp — 7 it) 7 n(AVjt) dx + / (a - ap) 7 n (AVu) dx 

J r 3 Jr 3 

= / 7 n (AVu)(TpXp - ap) dx 
J r 3 

= E 

& £ ^ext ? & CZ 1^3 

== E 

0"G^extj 0"CZr3 

We then deal edge by edge, considering two cases as in the proof of Theorem |2. 11 


/ (7 u{x a ) - a(x lT ))'y n (AVu) dx 

J U 

G a (u,v v ). 


In Case (i), where 7 u < a a.e. on a, we have G a (u,vv) = 0 since 7 n (AVx) = 0 in 
a. In Case (ii), we estimate G a (u, xp) using the Taylor expansion (5.12). □ 


6.2. The obstacle problem. With <jd £ L 2 (Q) an approximation of g. we con¬ 
sider the new convex set 


(6.3) 


/Cp = {x £ Aip o : U v v < g-D in 12} 


and we write the gradient scheme (2.7) with K. p instead of K, p. The following 


theorems are the equivalent for the obstacle problem of Theorems |6.1| and |6.2| 

Theorem 6.3 (Error estimates for the obsta cle problem with approximate bar¬ 
rier). Under the assumptions of Theorem\2.£\ if /Cp is not empty then there exists 


a unique solution u to the gradient scheme (12/7) in which /Cp has been replaced 


with /Cp. Moreover, Estimates (2.14) and (2.15) hold for any xp £ /Cp, provided 
that Ex>(u,vx>) is replaced with 

E v (u,vv) = E v (u,vv) + / (div(AVrt) + /)(<?p — g) dx. 

J n 

Proof. We follow exactly the proof of Theorem |2.8[ except that we introduce gp 


instead of g in (5.9). The first term in the right-hand side of this equation is then 


still bounded above by 0, and the second term is written 

[ (gv - nppp)(div(AVu) + /) dx = f (g v - g)(div(AVti) + /) dx 
Jo. Jn 

+ [ (g - IIpxp)(div(AVu) + /) dx. 

Jn 

The first term in this right-hand side corresponds to the additional term in Ep(u, xp), 
whereas the second term is the one handled in the proof of Theorem |2.8| □ 


Theorem 6.4 (Obstacle problem: order of convergence for piecewise constant 


reconstructions). LefD be a gradient discretisation, in the sense of Definition 2.5 
and let T = (A4,£,V) be a polytopal mesh of Cl. Let gx> £ T 2 (S2) be defined by 

\/K £ M , gv = 9 (xk) on K. 


Let u be the solution to (2.2), and let us assume that div(AVu) £ L 2 (UL) and that 
u — g £ W 2 ’ 00 ^). We also assume that there exists an interpolant xp £ Xp i0 such 
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that St> (u, vt>) < Cihj \4 , withCi not depending onV orT, and thatH'DVv = u(xk ) 
on K, for any K £ A4. 

Then K,x> 7 ^ 0 and, if u is the solution to \2.7\ with replaced with Kx>, 

(6.4) ||V-p7t — Vit|| i 2(Q)d + ||n-p7t — w||t2(q) < Ch.M + CWd(AVu) 

where C depends only on tt, A, u, g, C± and an upper bound ofCx>- 

Proof. The proof can be conducted by following the same ideas as in the proof of 
Theorem 16.21 □ 


Let us now compare our results with previous studies. In [7, 26], 0 (/i.m) error 
estimates are established for Pi and mixed finite elements applied to the Signorini 
problem and the obstacle problem. This order was obtained under the assumptions 
that u £ kF 1,00 (fi) D H 2 (Q) and a is constant for the Signorini problem and that 
A = Id, that u and g are in H 2 (Tl) for the case of the obstacle problem. Our results 
generalise these orders of convergence to the case of a Lipschitz-continuous A and 
a non-constant a (note that mixed finite elements are also part of the gradient 
schemes framework, see nans]). 

Studies of non-conforming methods for variational inequalities are more scarce. 
We cite [31], that applies Crouzeix-Raviart methods to the obstacle problem and 
obtains an order 0{Hm) under strong regularity assumptions, namely / £ L°°{Pl), 
u — g£ W 2,00 (f2), g £ H 2 (fl), A = Id and the free boundary has a finite length. For 
the Signorini problem, under the assumptions that a is constant, u £ VF 2,00 (f2) and 
A = I d , Wang and Hua [32] give a proof of an order 0(h M ) for Crouzeix-Raviart 
methods. The natural Crouzeix-Raviart interpolants are piecewise linear on the 
edges and the cells, and if we therefore use it as a-p for the Signorini problem, 
we see that under an H 2 regularity assumption on a, we have ||a — aT>\\L 2 (dn) < 
Hence, the additional term in Gti(u,vv) in Theorem 


Ch 2 M 


6.1 


is of order h 2 M and 


this theorem therefore gives back the known 0{Jim) order of convergence for the 
Crouzeix-Raviart scheme applied to the Signorini problem. This is obtained under 
slightly more general assumptions, since A does not need to be Id here. The same 
applies to the obstacle problem through Theorem |2.13| 

As shown in Section [3j our results also allowed us to obtain brand new orders of 
convergences, for methods not previously studied for (but relevant to) the Signorini 
or the obstacle problem. 

We finally notice that most of previous research investigates the Signorini prob¬ 
lem under the assumption that the barrier a is a constant. On the contrary, Theo¬ 
rems [2/7] [BTT] and [672] presented here are valid for non-constant barriers. 


7. Conclusion 

We developed a gradient scheme framework for two kinds of variational inequal¬ 
ities, the Signorini problem and the obstacle problem. We showed that several 
classical conforming or non-conforming numerical methods for diffusion equations 
fit into this framework, and we established general error estimates and orders of 
convergence for all methods in this framework, both for exact and approximate 
barriers. 

These results present new, and simpler, proofs of optimal orders of convergence 
for some methods previously studied for variational inequalities. They also allowed 
us to establish new orders of convergence for recent methods, designed to deal with 
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anisotropic heterogeneous diffusion PDEs on generic grids but not yet studied for 
the obstacle or the Signorini problem. 

We designed an HMM method for variational inequalities, and showed through 
numerical tests (including a new test case with analytical solution) the efficiency of 
these method. 

Acknowledgement: the authors thank Claire Chainais-Hillairet for providing 
the initial MATLAB HMM code. 


8. Appendix 

We recall here some classical notions on the normal trace of a function in H^ v . 
Let C be a Lipschtiz domain. Then there exists a surjective trace mapping 
7 : For if) € H^w, the normal trace 7 n(V’) £ (if 1 / 2 (df l))' of 

ip is defined by 

(8.1) ( 7 n (V>), 7 (u;)) = / ip-Vwdx + / dWip-wdx, for any w € 

Jo. Jq 

where (-, •) denotes the duality product between and (if ^(dfl))'. Stoke’s 

formula in if)] ( 0 ) x ifdiv(^) shows that this definition makes sense, i.e. the right- 
hand side is unchanged if we apply to w £ ii 1 (n) such that 7 (ie) = 7 (w). 
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